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Abstract -With Monte Carlo simulations, we study the creep motion of a domain wall in the two- 
dimensional random-field Ising model with a driving field. We observe the nonlinear field- velocity 
relation, and determine the creep exponent fi. To further investigate the universality class of the 
creep motion, we also measure the roughness exponent £ and energy barrier exponent ijj from the 
zero-field relaxation process. For strong disorder, the exponents are consistent with those of the 
Edwards- Wilkinson equation; for weak disorder, a different universality class is detected. 



Introduction. — In recent years the dynamics of elas- 
tic systems in disordered media has been a focus of the- 
oretical and experimental studies. Examples are charge- 
density waves, vortex lattices, domain walls in ferromag- 
netic or ferroelectric materials, contact lines, and fluid 
invasion in porous media [TUB]. In particular, the mag- 
netic domain- wall dynamics is an important topic in mag- 
netic devices, nanomaterials, thin films, and semiconduc- 
tors [71-110] . At zero temperature, the domain- wall motion 
exhibits a dcpinning transition. Due to the energy barri- 
ers created by the disorder, the domain wall is pinned and 
the velocity of the domain wall remains zero up to a criti- 
cal field H c [TTh14| . At finite temperature, the dcpinning 
transition is softened and the energy barriers can always 
be overcome by thermal activation. Thus, with a suffi- 
ciently small field, H <C H c , the domain- wall motion can 
reach a steady state, which is known as the creep motion 

Up to date, most theoretical approaches to the domain- 
wall dynamics in ferromagnetic materials are typically 
based on the Edwards- Wilkinson equation with quenched 
disorder (QEW). This equation is a phenomcnological 
model, and detailed microscopic structures and interac- 
tions of real materials are not concerned. To further un- 
derstand the domain-wall motion from a more fundamen- 
tal level, we should build lattice models based on micro- 
scopic structures and interactions. The random-field Ising 



model with a driving field (DRFIM) is a candidate, at 
least to capture robust features of the domain-wall mo- 
tion, although it does not include all interactions in real 
materials. It may go beyond the QEW equation, reveal 
new universality classes, explain a wider range of phys- 
ical phenomena, and allow a closer comparison with ex- 
periments. For example, the depinning transition of the 
two-dimensional DRFIM model has been examined with 
the short-time dynamic approach [THHH] • The critical ex- 
ponents are accurately determined, and the results indi- 
cate that the DRFIM model and QEW equation are not 
in a same universality class. In a very recent work, it is 
demonstrated that the relaxation state and relaxation-to- 
creep transition of the domain-wall motion can be simu- 
lated with the DRFIM model, and the Cole-Cole plot of 
the complex susceptibility obtained in the experiments is 
nicely reproduced [20 1 121]. 

The main purpose of this paper is to study the creep 
motion of the domain wall in the two-dimensional DRFIM 
model, in comparison with the QEW equation and exper- 
iments. In the creep regime of the domain-wall motion, 
the field- velocity relation is expected to be 



v~exp[-{H c /HyiT\, 



(1) 
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where [i is the creep exponent, usually used to classify the 
universality class of the materials in experiments [2ll5ll2"2"] . 
Based on the QEW equation, two universality classes have 
been identified for the field-induced domain- wall motion in 
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two dimensions [T7JI23] : A* = 1/4 belongs to the random- 
bond disorder, whereas = 1 to the random-field disorder, 
respectively. In the past years, most experiments confirm 
this result. However, recent progress in experiments sug- 
gests that the degree of surface roughness may alter the 
universality class Different values of /x between 1 

and 1/4 are also reported 25 , 26 . Due to the uncertainty 
of the domain-wall dimensionality, it leads to conflicts in 
identifying the universality classes of the materials [25]l2"7] . 
These experimental results challenge our theoretical un- 
derstanding. In the QEW equation, the degree of disorder 
is thought to be irrelevant to the universality class. 

On the other hand, based on the QEW equation, one 
may theoretically derive the scaling law [151128] 



2C + d - 2 
2-C ' 



(2) 



where C is the equilibrium roughness exponent. In the nu- 
merical study of the QEW equation, however, a violation 
of the scaling law in Eq. (|2|) is observed for low temper- 
atures or with strong disorder [15] , which implies a more 
complicated physical scenario. Moreover, it is argued in 
Ref . [29] that instead of Eq. ([2]) , the creep exponent should 
obey the scaling law 



2-C 



(3) 



where ip is the energy barrier exponent without a driving 
field [H1130]. The scaling laws in Eqs. © and © will 
be equivalent if rp is identical to the ground-state-energy 
fluctuation exponent (9, 



2( + d-2. 



(4) 



Theoretically, there are still some controversies on the 
value of V [I3[29l[31]. 

In order to fully understand the creep motion and its 
universality class, it is important to directly determine the 
exponents ip and f. This can be done by independently 
simulating the zero-field relaxation process of a domain 
wall [T5l[3T] ■ Since the driving field is absent, the cor- 
related domain-wall segment of a length £ overcomes the 
energy barrier t/(£) only by thermal activation. Accord- 
ing to the Arrhenius law, the time scale of this activation 
process is the order of exp[U(£)/T]. It is believed that 
the energy barrier may scale as U(£) ~ for a large 
£. Therefor, at long times, the correlation length should 
grows by a logarithmic law 



(5) 



As the domain wall propagates, it continues to roughen 
as an interface. The roughness function of the domain 
interface is believed to scale as 14. 18,32 34 



u 2 (t) ~ W, (6) 
where C is the equilibrium roughness exponent. 



Early numerical study of the DRFIM model also in- 
tended to tackle the creep motion [TB]. However, due 
to the limited range of the driving fields, the validity of 
Eq. ([T]) is not conclusively confirmed. Further understand- 
ing of the scaling laws in Eqs. ([3]) and Q is also lacking. 
In particular and more importantly, in this paper, we in- 
vestigate the possible dependence of the universality class 
on the strength of the disorder. 

The model. — The two-dimensional DRFIM model 
is defined by the Hamiltonian 



H = -JJ2 SiSj - J^{hi + H)Si, 



(7) 



<ij> 



where Si — ±1 is an Ising spin on the square lattice. The 
quenched random field hi is uniformly distributed within 
an interval [—A, A], and H is a homogeneous driving field. 
In this paper, we take the coupling constant J = 1. Since 
the two-dimensional random-field Ising model does not 
have a long-range ordered phase, one expects that domains 
grow spontaneously in the bulk at a finite temperature. 
This may lead to ambiguity in defining a single domain 
interface. If it is at low enough temperatures, however, 
the dynamic evolution of the bulk is negligible. In our 
study, the temperature is restricted to T max = 0.67. 

Our simulations are performed on a lattice, with a linear 
size 2L in the x direction and L in the y direction. An- 
tiperiodic and periodic boundary conditions are adopted 
in the x and y directions, respectively. To eliminate the 
pinning effect irrelevant for the quenched disorder, we ro- 
tate the square lattice such that the initial domain wall 
orients in the (11) direction of the square lattice, as shown 

in Refs. dungnHUHg. 

The initial state is a semiordered state with a perfect 
domain wall in the y direction. After preparing the initial 
state, we randomly select a spin, and flip it if the total en- 
ergy decreases after flipping. A Monte Carlo time step is 
defined by 2L 2 single-spin flips. Simulations are performed 
at two temperatures T = 0.33 and 0.67 with lattice size 
L = 256 up to t max = 10 6 . For the strength of random 
fields, two typical values of A are chosen, i.e., A = 0.5 
and 1.2. Simulations of different L are also performed 
to confirm that the finite-size effects are already negligi- 
bly small. For the simulations of the creep motion, the 
domain- wall velocity is averaged by at least 1000 samples. 
For the zero-field relaxation process, the total samples are 
5000. Errors are estimated by dividing the samples into 
three or four subgroups. 

To study the creep motion and zero-field relaxation of 
the domain interface, we first introduce a line magnetiza- 
tion 



/ \ 1 
m{y,t) = - 



' L 

E 

x=l 



S X y(t) 



(8) 



Here S xy (t) denotes a spin at site {x,y). The height func- 
tion of the interface is defined as 



%>*) = 2 Mf. *) + !]■ 



(9) 
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With the height function at hand, the average velocity of 
the interface can be calculated 

where < • • • > includes the statistical average and average 
over y. 

With the height function h(y, t), the roughness function 
of the interface is defined as 

u 2 (t) = {h{y,tf)-{h{y,t)) 2 . (11) 

A more informative quantity is the height correlation func- 
tion [33] . 

C(r,t) = (lh(y + r 7 t)-h(y,t)} 2 ). (12) 

It describes both the spatial correlation of the height func- 
tion in the y direction and the growth of the domain in- 
terface in the x direction. 

For a sufficiently large lattice, i.e., the correlation length 
C L, we should observe the standard scaling behavior 
for the roughness function in Eq. ©, and that for the 
height correlation function [T4l[T8l l32ti34] . 

C(r,t)=2u J 2 (t)f(r/0. (13) 

Monte Carlo simulations. — To investigate the 
characteristic relation between the driving field H and 
domain-wall velocity v, we measure the velocity for a given 
driving field when the steady state is reached. Theoret- 
ically the creep regime is expected to be far below the 
critical point H c . For a smaller driving field H, however, 
it takes a longer time to reach the steady state. In our sim- 
ulations, the longest time is 10 6 Monte Carlo time steps. 

In Fig. [TJ we plot the velocity v as a function of 
(Hc/H^/T. By adjusting the exponent /x, one may ob- 
serve a linear curve. Thus the nonlinear field-velocity re- 
lation in Eq. ([1]) is verified in the two-dimensional DR- 
FIM model. For A = 1.2, we obtain n = 1.02(5) and 
0.95(8) at T = 0.67 and 0.33 respectively, in good agree- 
ment with the theoretical prediction fj. = 1 based on the 
QEW equation. The small deviation from the exponen- 
tial law for large fields suggests that the system leaves the 
creep regime. For A = 0.5, however, the fitted exponent is 
(i = 0.63(5) and 0.59(4) at T = 0.67 and 0.33 respectively, 
significantly smaller than the theoretical value /j, = 1 of the 
QEW equation. To further confirm our results, we may 
plot lnf against H c /H on a log-log scale. From slopes 
of the curves, we measure the exponent ji. Or, we may 
directly fit the numerical data to the formula in Eq. (TT]) 
to extract fi. All these techniques yield similar values for 
the exponent [i. 

In Ref . [16] , it is stated that the numerical data could fit 
to the field- velocity relation in Eq. (fTJ) in the limit /j, 0. 
In the inset of Fig. Q] we plot v against ln(H c /H)/T on a 
semi-log scale. A tendency of a linear behavior is observed 
for large driving fields, supporting the results in Ref. [16] , 



Therefore we believe that the field regime considered in 
Ref. [TB] is too large to observe the creep motion. 

To verify the scaling laws in Eqs. ([3]) and ((4]), one needs 
to determine the exponents tJj and £ independently. Thus 
we simulate the zero-field relaxation process of the domain 
wall. In the absence of the driving field H, the growth law 
of the correlation length £(t) plays a crucial role. To ex- 
tract the correlation length £(t), we use a data collapse 
technique based on the scaling behavior of the height cor- 
relation function. According to Eq. (p~3|) . we fix t' = 10 6 , 
and rescale r of another t to £(i')/£(£)r and C(r,t) to 
w(t')/u)(t)C(r,t). If the ratio of £(*')/£(*) is properly cho- 
sen, C(r, t) at different time t collapse onto a single curve. 
This is shown in the inset of Fig. [2] In other words, one 
may estimate the ratio £(t')/£(t) from the data collapse of 
C(r,t). Thus, up to a constant £(io), the growth law 
is extracted. Alternatively, one may also determine 
by fitting C(r,t) to the scaling function [33] 

C(r, t) = A[tanh(r/B)} 2a , (14) 

where B is interpreted as the correlation length at t. Both 
methods yield the same results. 

In Fig. [2] is shown as a function of t on a log- 
log plot. A crossover behavior is observed. The corre- 
lation length seems to be described by a power law at 
early times, then obviously slows down after the transient 
regime, which implies a logarithmic growth. To extract 
the energy barrier exponent ip, we fit the correlation length 
to Eq. ([5]) in the time regime to < t < 10 6 . By increasing 
to until the fitted parameters become stabilized, we ob- 
serve the logarithmic growth law for almost three decades 
of time, shown by dash lines in Fig. @. The exponent 
tp w 1 for A = 1.2 is consistent with the theoretical predic- 
tion based on the QEW equation. In contrast, we obtain 
ip « 0.7 for A = 0.5. 

With the correlation length £ (t) at hand, one may mea- 
sure the roughness exponent £ from Eq. ([5]). In Fig. [3J 
w 2 (t) is plotted against £(t) on a log-log scale. A crossover 
behavior is observed. In the large-£(f) regime, the slopes 
of the curves yields C = 1.03(3) and 0.85(2) for 5 = 1.2 and 
0.5, respectively. Here we remind oursclf that £(f) grows 
very slowly by a logarithmic law, therefore it remains less 
than 10 2 even when f reaches 10 6 Monte Carlo time steps. 

In Table [1] all the measurements of the exponents n, ip, 
C, 0, and "0/(2 — C) are listed. For comparison, the QEW 
equation predicts that all these exponents are equal to 1 
[HUBS]- For strong disorder, i.e., A = 1.2, the exponents 
are consistent with those of the QEW equation. For weak 
disorder, i.e., A = 0.5, however, the exponents are ob- 
viously different, indicating a different universality class. 
Especially, the exponent \x as 0.6 is between the random- 
bold universality class /i = 1/4 and the random-field uni- 
versality class fi — 1 predicted by the QEW equation. This 
may give hints for understanding the experiments which 
yield a similar exponent (x [24H27] . On the other hand, 
although the exponents for the weak disorder depart from 
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those for the strong disorder, the scaling laws in Eqs. ^ 
and (@| always hold, i.e., ip « 9 and n ~ ip/(2 — Q within 
statistical errors. 

Why do the exponents depend on the strength of the 
disorder? At zero temperature, for A < 1, the depinning 
phase transition is of first order, and the transition field 
H c = A, since no overhangs can be created, while for 
A > 1, the depinning phase transition is of second order, 
and the transition field H c < A, since overhangs can be 
generated [TTl[T8l[T9j|37] . For the creep motion at nonzero 
temperatures, the driving field H is much smaller than A 
and H c , and the thermal activation plays an important 
role in driving the domain wall to propagate. For A < 1, 
however, overhangs can hardly be generated at low tem- 
peratures. In other words, a stronger disorder may induce 
overhangs, and overhangs may change the degree of the 
interface roughness, and alter the universality class of the 
domain interface. 

In the small-£(i) regime, the roughness exponent is ap- 
parently smaller than that in the large-£(i) regime, for ex- 
ample, ( = 0.65(3) for A = 0.5, which is close to ( = 2/3 of 
the random-bond disorder [3"51[3T)] . Interestingly, as shown 
in Fig. [3l the curves of different temperatures with a same 
strength of the disorder collapse onto a single curve, which 
implies that the temperature only affects the growth law of 
the correlation length This suggests that there may 
exist a temperature-independent characteristic length L c 
in the dynamic relaxation process. When the correlation 
length is smaller than L c , £(t) grows by a power law 
and the roughness exponent ( = 0.65(3), which may be 
comparable with the random-bond universality. When the 
correlation length is larger than L c , the growth law 
of (,(t) crosses over to a logarithmic one. 

Conclusion. — To summarize, we have simulated the 
creep motion of a domain wall in the two-dimensional 
random-field Ising model with a driving field. The non- 
linear field-velocity relation in Eq. ([1]) is confirmed, and 
the creep exponent /i is determined. Meanwhile, we also 
measure the roughness exponent £ and energy barrier ex- 
ponent ip from the zero-field relaxation process. For strong 
disorder, all the exponents are consistent with those pre- 
dicted by the Edwards- Wilkinson equation; for weak dis- 
order, however, the exponents are very much different, 
indicating a different universality class. Although the ex- 
ponents change with the strength of the disorder, the scal- 
ing laws in Eqs. ([3]) and Q, and therefore also in Eq. (|2|), 
always hold. 
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Fig. 1: The domain- wall velocity v is plotted as a function of 
(H c I H)^ /T for different disorders and temperatures on a semi- 
log scale. Dash lines show an exponential behavior. For A = 
f .2, /i = 1.02(5) and 0.95(8) at T = 0.67 and 0.33 respectively; 
for A = 0.5, (i = 0.63(5) and 0.59(4) at T = 0.67 and 0.33 
respectively. In the inset, v is plotted against \n(H c / H)/T for 
A = 1.2 and T = 0.67 on a semi-log scale. Errors of the data 
points are smaller than the symbol sizes. 



Fig. 2: The correlation length £(t) is displayed for different 
disorders and temperatures on a log-log scale. The curves of 
A = 1.2 are shifted for clarity. Dashed lines represent loga- 
rithmic fits to Eq. J5J). For A = 1.2, Tp = 1.05(7) and 0.98(9) 
at T = 0.67 and 0.33 respectively. For A = 0.5, V = 0.69(5) 
and 0.65(6) at T = 0.67 and 0.33 respectively. In the inset, 
the height correlation function C(r,t) of A = 1.2 and T = 0.67 
at different times are displayed with solid lines. Symbols show 
the data collapse according to Eq. (|13H 
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Fig. 3: The roughness function lo (t) is plotted against £(t) on 
a log-log scale. 





A = 1.2 


A = 0.5 


T 


0.67 0.33 


0.67 0.33 


A 4 


1.02(5) 0.95(8) 


0.63(5) 0.59(4) 




1.05(7) 0.98(9) 


0.69(5) 0.65(6) 


c 


1.03(3) 


0.85(2) 


8 


1.06(6) 


0.71(4) 


V7(2-C) 


1.08(8) 1.01(9) 


0.60(5) 0.57(6) 



Table 1: Exponents for the two-dimensional DRFIM model. 
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